## Courting the President replication
##
## contender & non-contender analysis 
## Estimand: ATE 
## Dependent variable: presIdeoVote 
##
## Reproduce: 
##      - Table A1 
##
## Kevin Quinn and Guoer Liu
## 11/19/2024

# setwd("") # set at your current working directory
set.seed(91800198)


source("func.R")
library(haven)
library(WeightIt)
library(Matching)
library(cem)
library(data.table)
library(RColorBrewer)
library(lattice)

##


## #############################################################
##
## Courting the President Replication
##
##      - Dependent variable: presIdeoVote
##      - Estimand: ATE
##      - Dataset: Contender 
##
##
## #############################################################




## Contender Judges

mydata <- read_dta("Contender.dta")
mydata <- as.data.frame(mydata)


mydata <- mydata[, c("presIdeoVote", "treatFinal0", "judgeJCS", "presDist",
                     "panelDistJCS", "circmed", "sctmed", "coarevtc",
                     "casepub", "judge")]
mydata <- na.omit(mydata)

cov.names.sub <- c("judgeJCS", "presDist", "panelDistJCS", "circmed",
                   "sctmed", "coarevtc", "casepub")

treat.form <- "treatFinal0 ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc + casepub"


reg.general.form <- "presIdeoVote ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc + casepub"




estimator <- NULL
estimate <- NULL
SE <- NULL
Eff.n <- NULL
Eff.n.ratio <- NULL
DFBETA.min <- NULL
DFBETA.max <- NULL
DFBETA.min.who <- NULL
DFBETA.max.who <- NULL
extrap.treat <- NULL
extrap.control <- NULL


## #############################################################
## general regression (two regression models)
estimator <- c(estimator, "reg.general")
mydata.1 <- mydata[mydata$treatFinal0==1,]
mydata.0 <- mydata[mydata$treatFinal0==0,]

lm.out.1 <- lm(reg.general.form, data=mydata.1)
lm.out.0 <- lm(reg.general.form, data=mydata.0)

m1 <- predict(lm.out.1, newdata=mydata)
m0 <- predict(lm.out.0, newdata=mydata)

estimate <- c(estimate, mean(m1 - m0))


## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    mydata.bs.1 <- mydata.bs[mydata.bs$treatFinal0==1,]
    mydata.bs.0 <- mydata.bs[mydata.bs$treatFinal0==0,]

    lm.out.1 <- lm(reg.general.form, data=mydata.bs.1)
    lm.out.0 <- lm(reg.general.form, data=mydata.bs.0)
    
    m1 <- predict(lm.out.1, newdata=mydata.bs)
    m0 <- predict(lm.out.0, newdata=mydata.bs)

    estim.bs[b] <- mean(m1 - m0)
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
X1 <- as.matrix(mydata[mydata$treatFinal0==1, cov.names.sub])
X1 <- cbind(1, X1)
X0 <- as.matrix(mydata[mydata$treatFinal0==0, cov.names.sub])
X0 <- cbind(1, X0)
n1 <- nrow(X1)
n0 <- nrow(X0)
n <- n1 + n0
treat.vec <- c(rep(1, n1), rep(0, n0))
y1 <- mydata[mydata$treatFinal0==1, "presIdeoVote"]
y0 <- mydata[mydata$treatFinal0==0, "presIdeoVote"]
X <- rbind(X1, X0)
y <- as.matrix(c(y1, y0), nrow=n, ncol=1)
ones.n <- matrix(1, nrow=n, ncol=1)

w1 <- t(ones.n) %*% X %*% solve(t(X1) %*% X1) %*% t(X1)
w0 <- t(ones.n) %*% X %*% solve(t(X0) %*% X0) %*% t(X0)
w <- as.vector(c(w1, w0))

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treat.vec, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

mydata.reg.ate <- rbind(mydata.1, mydata.0)
mydata.reg.ate$w <- w
mydata.reg.ate$DFBETA.out <- DFBETA.out

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata.reg.ate[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata.reg.ate[max.ind, "judge"])


extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)



## END general regression (two regression models)
## #############################################################





## #############################################################
## general regression (constant effect)
estimator <- c(estimator, "reg.constant")
reg.cons.form <- "presIdeoVote ~ treatFinal0 + judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc + casepub"
lm.out.cons <- lm(reg.cons.form, data=mydata)

estimate <- c(estimate, coef(lm.out.cons)["treatFinal0"])

## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]
    lm.out.bs <- lm(reg.cons.form, data=mydata.bs)
    estim.bs[b] <- coef(lm.out.bs)["treatFinal0"]
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
X1 <- as.matrix(mydata[mydata$treatFinal0==1, cov.names.sub])
X1 <- cbind(1, X1)
X0 <- as.matrix(mydata[mydata$treatFinal0==0, cov.names.sub])
X0 <- cbind(1, X0)
n1 <- nrow(X1)
n0 <- nrow(X0)
n <- n1 + n0
treat.vec <- c(rep(1, n1), rep(0, n0))
y1 <- mydata[mydata$treatFinal0==1, "presIdeoVote"]
y0 <- mydata[mydata$treatFinal0==0, "presIdeoVote"]
X <- rbind(X1, X0)
y <- as.matrix(c(y1, y0), nrow=n, ncol=1)
U <- cbind(X, treat.vec)
U1 <- cbind(X, 1)
U0 <- cbind(X, 0)
ones.n <- matrix(1, nrow=n, ncol=1)

w1 <- t(ones.n) %*% U1 %*% solve(t(U) %*% U) %*% t(U) -
    t(ones.n) %*% U0 %*% solve(t(U) %*% U) %*% t(U)
w0 <- t(ones.n) %*% U0 %*% solve(t(U) %*% U) %*% t(U) -
    t(ones.n) %*% U1 %*% solve(t(U) %*% U) %*% t(U)

w <- w0
w[1:n1] <- w1[1:n1]
w <- as.vector(w)
w1 <- w[1:n1]
w0 <- w[(n1+1):n]

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treat.vec, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

mydata.reg_cons.ate <- rbind(mydata.1, mydata.0)
mydata.reg_cons.ate$w <- w
mydata.reg_cons.ate$DFBETA.out <- DFBETA.out

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata.reg_cons.ate[min.ind, "judge"]) 
DFBETA.max.who <- c(DFBETA.max.who, mydata.reg_cons.ate[max.ind, "judge"]) 

extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)



## END general regression (constant effect)
## #############################################################





## #############################################################
## propensity score weighting (logit)
estimator <- c(estimator, "sipw.logit")

glm.out <- glm(treat.form, data=mydata, family=binomial(logit))
mydata$pscores <- glm.out$fitted
mydata$w <- 1 / mydata$pscores
mydata$w[mydata$treatFinal0==0] <- 1 /
    (1 - mydata$pscores[mydata$treatFinal0==0])

ATE.hat <- sum(mydata$w * mydata$treatFinal0 * mydata$presIdeoVote) /
    sum(mydata$w * mydata$treatFinal0) -
    sum(mydata$w * (1-mydata$treatFinal0) * mydata$presIdeoVote) /
    sum(mydata$w * (1-mydata$treatFinal0))

estimate <- c(estimate, ATE.hat)


## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    glm.out.bs <- glm(treat.form, data=mydata.bs,
                      family=binomial(logit))
    mydata.bs$pscores <- glm.out.bs$fitted
    mydata.bs$w <- 1 / mydata.bs$pscores
    mydata.bs$w[mydata.bs$treatFinal0==0] <- 1 /
        (1 - mydata.bs$pscores[mydata.bs$treatFinal0==0])
    
    ATE.hat.bs <- sum(mydata.bs$w * mydata.bs$treatFinal0 * mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * mydata.bs$treatFinal0) -
        sum(mydata.bs$w * (1-mydata.bs$treatFinal0) * mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * (1-mydata.bs$treatFinal0))
    
    
    estim.bs[b] <- ATE.hat.bs
}
SE <- c(SE, sd(estim.bs))


w0 <- mydata$w[mydata$treatFinal0==0]
w1 <- mydata$w[mydata$treatFinal0==1]
n1 <- sum(mydata$treatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$treatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)


mydata.sipw.logit.ate <- mydata
mydata.sipw.logit.ate$DFBETA.out <- DFBETA.out

## END propensity score weighting (logit)
## #############################################################









## #############################################################
## entropy balancing
estimator <- c(estimator, "entropy.balancing")


eb.out <- weightit(as.formula(treat.form), data=mydata, method="ebal",
                   estimand="ATE")
mydata$w <- eb.out$weights

wls.out <- lm(presIdeoVote ~ treatFinal0, weights=w, data=mydata)

ATE.hat <- coef(wls.out)[2]

estimate <- c(estimate, ATE.hat)
SE <- c(SE, sqrt(vcov(wls.out)[2,2]))


w0 <- mydata$w[mydata$treatFinal0==0]
w1 <- mydata$w[mydata$treatFinal0==1]
n1 <- sum(mydata$treatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$treatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)


mydata.entropy.balancing.ate <- mydata
mydata.entropy.balancing.ate$DFBETA.out <- DFBETA.out



## END entropy balancing
## #############################################################





## #############################################################
## genetic matching
estimator <- c(estimator, "NN.matching")

treat.vec <- mydata$treatFinal0
X <- as.matrix(mydata[, c(cov.names.sub)])
y <- mydata$presIdeoVote

# Matching by linearized propensity score

glm.out <- glm(treat.vec ~ X, family=binomial(logit))
l.pscores <- predict(glm.out)


match.out <- Match(Y = y, Tr=treat.vec, X=l.pscores, estimand="ATE",
                   M=1, ties=FALSE, replace=TRUE, BiasAdjust=FALSE)


estimate <- c(estimate, match.out$est)
SE <- c(SE, match.out$se.standard)

mydata$w <- 0
for (i in 1:nrow(mydata)){
     mydata$w[i] <- sum(match.out$index.control == i) +
            sum(match.out$index.treated == i)
}

w0 <- mydata$w[mydata$treatFinal0==0]
w1 <- mydata$w[mydata$treatFinal0==1]
n1 <- sum(mydata$treatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$treatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)

mydata.NN.matching.ate <- mydata


## END genetic matching
## #############################################################


con.out.data.ate <- data.frame(estimator=estimator,
                       estimate=estimate,
                       SE=SE,
                       #Z = estimate / SE,
                       Eff.n=Eff.n,
                       Eff.n.ratio=Eff.n.ratio,
                       DFBETA.min=DFBETA.min,
                       DFBETA.min.who=DFBETA.min.who,
                       DFBETA.max=DFBETA.max,
                       DFBETA.max.who=DFBETA.max.who,
                       extrap.control=extrap.control,
                       extrap.treat=extrap.treat)

print(con.out.data.ate)


#library(xtable)
#xtable(con.out.data.ate, digits=c(0, 0, rep(3, 2), 1, rep(3, 2), 0, 3, 0, rep(3, 2)))




## #############################################################
##
## Courting the President Replication
##
##      - Dependent variable: presIdeoVote
##      - Estimand: ATE
##      - Dataset: Non-Contender 
##
##
## #############################################################


set.seed(91800198)


## Non-Contender Judges

mydata <- read_dta("NonContender.dta")
mydata <- as.data.frame(mydata)

## drop casepub b/c no variation in noncontender data
mydata <- mydata[, c("presIdeoVote", "nSLtreatFinal0",
                     "judgeJCS", "presDist",
                     "panelDistJCS", "circmed", "sctmed", "coarevtc",
                     "judge")]
mydata <- na.omit(mydata)

cov.names.sub <- c("judgeJCS", "presDist", "panelDistJCS", "circmed",
                   "sctmed", "coarevtc")

treat.form <- "nSLtreatFinal0 ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"


reg.general.form <- "presIdeoVote ~ judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"




estimator <- NULL
estimate <- NULL
SE <- NULL
Eff.n <- NULL
Eff.n.ratio <- NULL
DFBETA.min <- NULL
DFBETA.max <- NULL
DFBETA.min.who <- NULL
DFBETA.max.who <- NULL
extrap.treat <- NULL
extrap.control <- NULL


## #############################################################
## general regression (two regression models)
estimator <- c(estimator, "reg.general")
mydata.1 <- mydata[mydata$nSLtreatFinal0==1,]
mydata.0 <- mydata[mydata$nSLtreatFinal0==0,]

lm.out.1 <- lm(reg.general.form, data=mydata.1)
lm.out.0 <- lm(reg.general.form, data=mydata.0)

m1 <- predict(lm.out.1, newdata=mydata)
m0 <- predict(lm.out.0, newdata=mydata)

estimate <- c(estimate, mean(m1 - m0))


## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    mydata.bs.1 <- mydata.bs[mydata.bs$nSLtreatFinal0==1,]
    mydata.bs.0 <- mydata.bs[mydata.bs$nSLtreatFinal0==0,]

    lm.out.1 <- lm(reg.general.form, data=mydata.bs.1)
    lm.out.0 <- lm(reg.general.form, data=mydata.bs.0)
    
    m1 <- predict(lm.out.1, newdata=mydata.bs)
    m0 <- predict(lm.out.0, newdata=mydata.bs)

    estim.bs[b] <- mean(m1 - m0)
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
X1 <- as.matrix(mydata[mydata$nSLtreatFinal0==1, cov.names.sub])
X1 <- cbind(1, X1)
X0 <- as.matrix(mydata[mydata$nSLtreatFinal0==0, cov.names.sub])
X0 <- cbind(1, X0)
n1 <- nrow(X1)
n0 <- nrow(X0)
n <- n1 + n0
treat.vec <- c(rep(1, n1), rep(0, n0))
y1 <- mydata[mydata$nSLtreatFinal0==1, "presIdeoVote"]
y0 <- mydata[mydata$nSLtreatFinal0==0, "presIdeoVote"]
X <- rbind(X1, X0)
y <- as.matrix(c(y1, y0), nrow=n, ncol=1)
ones.n <- matrix(1, nrow=n, ncol=1)

w1 <- t(ones.n) %*% X %*% solve(t(X1) %*% X1) %*% t(X1)
w0 <- t(ones.n) %*% X %*% solve(t(X0) %*% X0) %*% t(X0)
w <- as.vector(c(w1, w0))

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treat.vec, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

mydata.noncontender.reg.ate <- rbind(mydata.1, mydata.0)
mydata.noncontender.reg.ate$w <- w
mydata.noncontender.reg.ate$DFBETA.out <- DFBETA.out

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, 
                    mydata.noncontender.reg.ate[min.ind, "judge"])

DFBETA.max.who <- c(DFBETA.max.who, 
                    mydata.noncontender.reg.ate[max.ind, "judge"])


extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)



## END general regression (two regression models)
## #############################################################





## #############################################################
## general regression (constant effect)
estimator <- c(estimator, "reg.constant")
reg.cons.form <- "presIdeoVote ~ nSLtreatFinal0 + judgeJCS + presDist + panelDistJCS + circmed + sctmed + coarevtc"
lm.out.cons <- lm(reg.cons.form, data=mydata)

estimate <- c(estimate, coef(lm.out.cons)["nSLtreatFinal0"])

## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]
    lm.out.bs <- lm(reg.cons.form, data=mydata.bs)
    estim.bs[b] <- coef(lm.out.bs)["nSLtreatFinal0"]
}
SE <- c(SE, sd(estim.bs))


## get the WATE weights
X1 <- as.matrix(mydata[mydata$nSLtreatFinal0==1, cov.names.sub])
X1 <- cbind(1, X1)
X0 <- as.matrix(mydata[mydata$nSLtreatFinal0==0, cov.names.sub])
X0 <- cbind(1, X0)
n1 <- nrow(X1)
n0 <- nrow(X0)
n <- n1 + n0
treat.vec <- c(rep(1, n1), rep(0, n0))
y1 <- mydata[mydata$nSLtreatFinal0==1, "presIdeoVote"]
y0 <- mydata[mydata$nSLtreatFinal0==0, "presIdeoVote"]
X <- rbind(X1, X0)
y <- as.matrix(c(y1, y0), nrow=n, ncol=1)
U <- cbind(X, treat.vec)
U1 <- cbind(X, 1)
U0 <- cbind(X, 0)
ones.n <- matrix(1, nrow=n, ncol=1)

w1 <- t(ones.n) %*% U1 %*% solve(t(U) %*% U) %*% t(U) -
    t(ones.n) %*% U0 %*% solve(t(U) %*% U) %*% t(U)
w0 <- t(ones.n) %*% U0 %*% solve(t(U) %*% U) %*% t(U) -
    t(ones.n) %*% U1 %*% solve(t(U) %*% U) %*% t(U)

w <- w0
w[1:n1] <- w1[1:n1]
w <- as.vector(w)
w1 <- w[1:n1]
w0 <- w[(n1+1):n]

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(y, treat.vec, w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

mydata.noncontender.reg_cons.ate <- rbind(mydata.1, mydata.0)
mydata.noncontender.reg_cons.ate$w <- w
mydata.noncontender.reg_cons.ate$DFBETA.out <- DFBETA.out

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, 
                    mydata.noncontender.reg_cons.ate[min.ind, "judge"]) 
DFBETA.max.who <- c(DFBETA.max.who, 
                    mydata.noncontender.reg_cons.ate[max.ind, "judge"]) 

extrap.1 <- sum(abs(w1) * (w1 < 0)) / sum(w1 * (w1 >= 0))
extrap.0 <- sum(abs(w0) * (w0 < 0)) / sum(w0 * (w0 >= 0))

extrap.treat <- c(extrap.treat, extrap.1)
extrap.control <- c(extrap.control, extrap.0)


## END general regression (constant effect)
## #############################################################





## #############################################################
## propensity score weighting (logit)
estimator <- c(estimator, "sipw.logit")

glm.out <- glm(treat.form, data=mydata, family=binomial(logit))
mydata$pscores <- glm.out$fitted
mydata$w <- 1 / mydata$pscores
mydata$w[mydata$nSLtreatFinal0==0] <- 1 /
    (1 - mydata$pscores[mydata$nSLtreatFinal0==0])

ATE.hat <- sum(mydata$w * mydata$nSLtreatFinal0 * mydata$presIdeoVote) /
    sum(mydata$w * mydata$nSLtreatFinal0) -
    sum(mydata$w * (1-mydata$nSLtreatFinal0) * mydata$presIdeoVote) /
    sum(mydata$w * (1-mydata$nSLtreatFinal0))

estimate <- c(estimate, ATE.hat)


## SE via nonparametric bootstrap
B <- 5000
n <- nrow(mydata)
estim.bs <- rep(NA, B)
for (b in 1:B){
    bs.inds <- sample(1:n, n, replace=TRUE)
    mydata.bs <- mydata[bs.inds,]

    glm.out.bs <- glm(treat.form, data=mydata.bs,
                      family=binomial(logit))
    mydata.bs$pscores <- glm.out.bs$fitted
    mydata.bs$w <- 1 / mydata.bs$pscores
    mydata.bs$w[mydata.bs$nSLtreatFinal0==0] <- 1 /
        (1 - mydata.bs$pscores[mydata.bs$nSLtreatFinal0==0])
    
    ATE.hat.bs <- sum(mydata.bs$w * mydata.bs$nSLtreatFinal0 * mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * mydata.bs$nSLtreatFinal0) -
        sum(mydata.bs$w * (1-mydata.bs$nSLtreatFinal0) * mydata.bs$presIdeoVote) /
        sum(mydata.bs$w * (1-mydata.bs$nSLtreatFinal0))
    
    
    estim.bs[b] <- ATE.hat.bs
}
SE <- c(SE, sd(estim.bs))


w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)


mydata.noncontender.sipw.logit.ate <- mydata
mydata.noncontender.sipw.logit.ate$DFBETA.out <- DFBETA.out

## END propensity score weighting (logit)
## #############################################################









## #############################################################
## entropy balancing
estimator <- c(estimator, "entropy.balancing")


eb.out <- weightit(as.formula(treat.form), data=mydata, method="ebal",
                   estimand="ATE")
mydata$w <- eb.out$weights

wls.out <- lm(presIdeoVote ~ nSLtreatFinal0, weights=w, data=mydata)

ATE.hat <- coef(wls.out)[2]

estimate <- c(estimate, ATE.hat)
SE <- c(SE, sqrt(vcov(wls.out)[2,2]))


w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)


mydata.noncontender.entropy.balancing.ate <- mydata
mydata.noncontender.entropy.balancing.ate$DFBETA.out <- DFBETA.out



## END entropy balancing
## #############################################################





## #############################################################
## genetic matching
estimator <- c(estimator, "NN.matching")

treat.vec <- mydata$nSLtreatFinal0
X <- as.matrix(mydata[, c(cov.names.sub)])
y <- mydata$presIdeoVote

# Matching by linearized propensity score

glm.out <- glm(treat.vec ~ X, family=binomial(logit))
l.pscores <- predict(glm.out)


match.out <- Match(Y = y, Tr=treat.vec, X=l.pscores, estimand="ATE",
                   M=1, ties=FALSE, replace=TRUE, BiasAdjust=FALSE)


estimate <- c(estimate, match.out$est)
SE <- c(SE, match.out$se.standard)

mydata$w <- 0
for (i in 1:nrow(mydata)){
     mydata$w[i] <- sum(match.out$index.control == i) +
            sum(match.out$index.treated == i)
}

w0 <- mydata$w[mydata$nSLtreatFinal0==0]
w1 <- mydata$w[mydata$nSLtreatFinal0==1]
n1 <- sum(mydata$nSLtreatFinal0)
n <- nrow(mydata)
n0 <- n - n1

n.eff.0 <- (sum(w0)^2) / sum(w0^2)
n.eff.1 <- (sum(w1)^2) / sum(w1^2)
n.eff <- (n.eff.1 * n.eff.0) / (n.eff.1 + n.eff.0)
n.eff.max <- (n1 * n0) / (n1 + n0)
Eff.n <- c(Eff.n, n.eff)
Eff.n.ratio <- c(Eff.n.ratio, n.eff/n.eff.max)

DFBETA.out <- DFBETA.weighted(mydata$presIdeoVote, mydata$nSLtreatFinal0,
                              mydata$w)
DFBETA.min <- c(DFBETA.min, min(DFBETA.out))
DFBETA.max <- c(DFBETA.max, max(DFBETA.out))

min.ind <- which.min(DFBETA.out)
max.ind <- which.max(DFBETA.out)

DFBETA.min.who <- c(DFBETA.min.who, mydata[min.ind, "judge"])
DFBETA.max.who <- c(DFBETA.max.who, mydata[max.ind, "judge"])


extrap.treat <- c(extrap.treat, NA)
extrap.control <- c(extrap.control, NA)

mydata.noncontender.NN.matching.ate <- mydata


## END genetic matching
## #############################################################


noncon.out.data.ate <- data.frame(estimator=estimator,
                       estimate=estimate,
                       SE=SE,
                       #Z = estimate / SE,
                       Eff.n=Eff.n,
                       Eff.n.ratio=Eff.n.ratio,
                       DFBETA.min=DFBETA.min,
                       DFBETA.min.who=DFBETA.min.who,
                       DFBETA.max=DFBETA.max,
                       DFBETA.max.who=DFBETA.max.who,
                       extrap.control=extrap.control,
                       extrap.treat=extrap.treat)

print(noncon.out.data.ate)


#library(xtable)
#xtable(noncon.out.data.ate, digits=c(0, 0, rep(3, 5), 0, 3, 0, rep(3, 2)))

